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DOUBLE-BASE  PROPELLANTS 


Stephan  R.  Bilyk  and  Michael  J.  Scheidler 


Weapons  and  Materials  Research  Directorate 
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Abstract.  The  intense  shearing  that  occurs  in  propellants  during  impulsive  loading  can  lead 
to  initiation.  In  an  effort  to  determine  useful  shear  initiation  criteria,  the  U.S.  Army 
Research  Laboratory  has  developed  a  dynamic  shear  punch  test  using  a  modified  split- 
Hopkinson  bar.  Varying  the  striker  bar’s  velocity  and  length  controls  the  shear  rate  and 
duration.  Shear  velocities  approaching  100  m/s  and  durations  as  long  as  0.2  ms  are  possible. 
Experimental  results  have  been  obtained  for  several  energetic  materials  and  a  nonreacting 
polymer,  polycarbonate  (PC).  This  paper  presents  a  detailed  analysis  used  to  obtain 
constitutive  behavior  and  shear  initiation  for  double-base  propellants  and  computational 
results  of  the  shear  punch  test.  For  the  simulations,  the  viscoSCRAM  constitutive  model 
was  used  to  describe  viscoelasticity,  cracking  and  ignition  in  the  propellant  when  subjected 
to  dynamic  shear  loading  conditions.  First,  we  will  present  the  analysis  used  to  obtain 
viscoelastic  material  parameters.  The  stress  relaxation  function  for  the  linear  viscoelastic 
response  was  obtained  by  using  time-temperature  superposition  to  generate  a  master  curve 
from  Dynamic  Mechanical  Analysis  (DMA)  data.  Next,  the  effect  of  initial  crack  size  and 
critical  hot  spot  duration  on  the  ignition  threshold  temperature  was  examined.  The  validity 
of  the  constitutive  relation,  failure  criterion,  and  shear  initiation  is  determined  based  on 
their  ability  to  predict  the  observed  response  from  the  dynamic  shear  punch  test. 


INTRODUCTION 

Energetic  materials  are  often  ranked  in  terms  of 
their  sensitivity  when  subjected  to  shock,  shear,  and 
thermal  stimuli.  The  goal  for  military  applications 
is  to  develop  initiation  criteria  under  each  stimuli  as 
well  as  a  fundamental  understanding  of  coupled 
behavior.  Several  useful  analytical  models  and 
experiments  already  exist  for  shock  and  thermal 
stimuli.  However,  initiation  due  to  shear  loading  is 
complex  and  poorly  understood.  Many  hazardous 
scenarios  such  as  hot  metal  fragments  impacting  an 
explosive  canister  can  lead  to  shear  initiation  of  an 
energetic.  Shear  initiation  occurs  at  timescales  over 


tens  or  hundreds  of  microseconds,  an  order  of 
magnitude  larger  than  shock  loading.  Energy  is 
deposited  in  localized  regions  causing  a  local 
temperature  rise,  which  for  some  energetics  can 
even  lead  to  the  development  of  adiabatic  shear 
bands. 

It  is  generally  accepted  that  initiation  of  an 
energetic  is  a  thermal  process1.  High  pressure 
accelerates  chemical  reactions,  but  most  often  does 
not  initiate  them.  Therefore,  critical  factors  to 
initiate  reaction  are  those  that  generate  heat  by 
direct  application  or  by  the  conversion  of 
mechanical  or  electrical  energy  to  heat.  This  paper 
describes  non-shock  mechanical  stimuli  sufficient 


to  create  local  regions,  so  called  “hot  spots”,  which 
can  lead  to  thermal  ignition.  By  non-shock  ignition 
we  mean  that  there  is  an  energy  release  but  no 
shock  wave. 

For  thermal  ignition  due  to  mechanical 
stimulus,  it  is  not  necessary  to  heat  the  bulk  of  the 
energetic  since  the  locally  created  hot  spots  may 
reach  sufficiently  high  temperatures.  Energetic 
materials  are  a  heterogeneous  mixture  of 
polycrystalline  explosive,  binder,  and  additives 
including  voids  created  during  material  processing. 
Mechanical  loading  can  nucleate  hot  spots 
(commonly  in  void  regions)  but  only  a  few  become 
critical  hot  spots.  These  critical  hot  spots  ignite  the 
energetic  if  the  generation  of  heat  in  the  localized 
volume  is  greater  than  the  heat  lost  to  the 
surroundings.  In  their  monograph  research  on  the 
topic,  Bowden  and  Joffe2  estimated  critical  hot  spot 
parameters  as  typically  of  micron  size  (0.1  -10pm), 
lasting  for  lOps  to  1ms,  and  reaching  temperatures 
of  approximately  700K.  Clearly,  if  local 
temperatures  are  high,  the  size  and  duration  can  be 
smaller.  Hot  spots  form  during  the  interaction  of 
stress  waves  with  material  defects  and  depend  on 
the  mechanical,  thermal  and  chemical  properties  of 
the  energetic.  There  are  different  mechanisms  at  the 
microstructural  length  scale  that  can  create  hot  spot 
ignition.  These  include  jetting  of  material  grains, 
hydrodynamic  pore  collapse,  viscous  heating,  shear 
localization,  friction  between  grains,  internal  shear 
and  shock  interaction  with  second  phase  particles3,4. 
The  dominant  mechanism  for  producing  the  hot 
spot  has  not  been  generally  agreed  upon.  However, 
Dienes5  analytically  showed  that  the  largest 
contribution  to  potential  heat  generation  is  the 
frictional  forces  on  shear  crack  surfaces. 

Double  base  propellants  are  composed  of 
nitrocellulose  and  stirred  with  a  reactive  plasticizer 
liquid  nitrate  ester  such  as  nitroglycerine  which  also 
affects  the  oxygen  balance.  Stabilizers  and 
gelatinizers  are  often  added  and  the  paste  is  hot 
rolled  processed  and  pressed  without  the  use  of  a 
solvent.  The  plasticizer  is  used  to  adjust  the  oxygen 
balance  which  affects  the  energy  output  and 
reaction  temperature6.  This  class  of  propellant 
powders  is  often  used  in  large  caliber  guns  and 
solid  rockets. 

Initially,  the  activator  punch  test  was 
developed  to  study  shear  initiation7.  This  test  was 
limited  since  it  was  difficult  to  control  the  shear 
velocity  independently  of  the  pressure  and  the 


pressure  on  the  shear  surface  was  not  well  known. 
Recently,  Krzewinski  et  al.8  developed  a  shear 
punch  test  at  ARL.  The  shear  punch  test  uses  a 
modified  Kolsky  bar  technique  and  obtains  data  for 
shear  initiation  of  energetic  materials  subjected  to 
dynamic  loading  conditions.  In  addition,  some  non- 
energetic  polymer  materials  such  as  polycarbonate 
(PC)  have  been  used  as  surrogate  specimens  for 
comparison  purposes. 

This  paper  first  establishes  the  shear  punch 
test,  describes  the  constitutive  model  used  for  the 
propellant  and  finally  provides  computational 
results  for  the  shear  punch  as  validation  of  the 
energetic  material  response  model. 

EXPERIMENTAL  DESCRIPTION 

The  apparatus  used  for  the  shear  punch 
test  developed  by  Krzewinski  et  al.8  was  a  modified 
Kolsky  bar,  as  shown  in  Figure  1.  The  striker, 
incident,  and  output  bars  were  1.27cm  diameter 
350-maraging  steel.  The  incident  and  output  bars 
were  150cm  in  length,  while  the  striker  bar  was 
available  in  25,  50,  and  55cm  lengths.  The  varying 
striker  lengths  gave  nominal  pulse  durations  of  100, 
200,  and  220ps,  respectively8.  The  specimen  had  a 
diameter  of  1 .905cm  and  a  length  of  1 .27cm.  As  the 
compression  wave  travels  during  the  entire  test,  the 
striker,  input  and  output  bars  and  the  holder  remain 
elastic.  The  specimen  is  the  only  material  that 
undergoes  plastic  deformation  but  not  at  constant 
strain  rate. 

The  experimental  measurements  are  also 
shown  (boxed)  in  Figure  1.  Impact  velocity  was 
measured  using  three  fiber  optic  wires  and  an 
optical  detector.  Two  strain  gages  were  mounted 
near  the  center  of  the  input  and  output  bars  to 
measure  the  incident,  reflected,  and  transmitted 
strains.  Finally,  a  scanning  electron  microscope 
(SEM)  was  used  to  measure  the  punch  and  dent 
displacements  of  the  specimen  as  well  as  examine 
any  fracture  regions. 


FIGURE  1.  Schematic  of  Shear  Punch  Test  and 
data  collection  (Not  to  Scale). 


A  special  shock  absorber  and  transfer 
piston  (not  shown  in  Figure  1)  were  designed  to 
prevent  reverse  bar  motion  whenever  the  specimen 
reacted  violently.  Thin  polyethylene  disks  were  also 
placed  between  the  specimen  and  incident/output 
bars  for  impedance  matching.  Copper  (3mil)  and 
Kaptan  (5mil)  disks  were  placed  between  the  striker 
and  incident  bar  to  reduce  ringing  and  wave  shape  a 
nearly  rectangular  incident  compressive  pulse.  The 
specimen  holder  was  made  from  17-4  PH  stainless 
steel  and  consisted  of  three  pieces  held  together 
with  six  high-strength  bolts.  In  addition,  vacuum 
grease  was  applied  between  the  specimen  and 
specimen  holder  to  fill  any  voids  and  reduce 
friction  at  the  interfaces.  With  the  applied  grease, 
one  can  conclude  that  all  initiations  occurred 
because  of  the  shearing  within  the  specimen. 

A  typical  deformed  specimen  shape 
observed  by  Krzewinski  et  al.8  is  shown  in  Figure 
2.  The  specimen  shown  is  a  double-base  propellant, 
PI.  Note  also  in  Figure  2,  that  the  shear  surface  has 
localized  and  runs  along  the  outer  radial  edges  of 
the  incident  bar.  For  this  dynamic  test,  the  loading 
on  the  P 1  specimen  was  great  enough  to  eventually 
fracture  the  specimen  along  the  shear  surface. 


^  —  Punch 

t 
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FIGURE  2.  Typical  specimen  deformation  and 
idealized  shear  surface  (dotted  line)8. 

DESCRIPTION  OF  THE  CONSTITUTIVE 
MODEL 

For  the  simulations,  the  constitutive 
behavior  of  the  specimen  was  modeled  using 
viscoSCRAM9.  The  model  captures  rate 
dependence  (linear- viscoelastic),  damage 
accumulation  (statistical-crack-mechanics), 

adiabatic  mechanical  heating  and  chemical  heating 
that  are  apparent  for  some  energetics.  Furthermore, 
the  model  does  not  include  heat  conduction  because 
it  is  too  slow  compared  to  the  deformation  time 
scale. 

The  mechanical  response  has  two 
constitutive  assumptions.  The  first  is  that  the  strain 
rate  can  be  decoupled  into  viscoelastic  and 
deviatoric  material  damage  components.  The 
second  is  that  the  shear  stress  is  determined  from 


the  viscoelastic  strain  rate.  The  viscoelastic  portion 
is  based  on  the  work  of  Addessio  and  Johnson10  and 
the  damage  model  uses  the  statistical  crack 
mechanics  (SCRAM)  approach  of  Dienes11.  The 
ignition  criterion  is  based  on  whether  the  localized 
temperature  increase  due  to  friction  between  sliding 
crack  faces  is  sufficient  for  the  hot  spot  to  become 
critical.  The  following  sub-sections  describe  each 
sub-model  and  relevant  material  parameters 
required  for  each  sub-model. 

Viscoelastic  Parameters 

For  sufficiently  small  strains,  the 

propellant  may  be  modeled  as  a  linear  viscoelastic 
solid.  Assuming  that  the  volumetric  response  is 
purely  elastic,  the  viscoelastic  response  is 

characterized  by  the  stress  relaxation  function  in 
shear,  denoted  by  G.  This  function  can  be 
approximated  by  a  sum  of  decaying  exponentials, 
referred  to  as  a  Prony  series: 

G(t)  =  G„+fjGne-th»  (1) 

n= 1 

Here  t  denotes  time,  xn  are  the  relaxation  times,  and 
Goo=  G(oo)  is  the  equilibrium  elastic  modulus.  The 
coefficients  Gn  are  also  moduli,  sometimes  referred 
to  as  spring  constants  since  this  relation  may  be 
obtained  by  analogy  with  spring  and  dashpot 

models.  The  relation  (1)  is  also  known  as  a 

generalized  Maxwell  model.  The  instantaneous 
elastic  modulus  is  the  value  of  the  stress  relaxation 
function  at  time  zero: 

G(0)  =  Goo+fjGn  (2) 

n= 1 

It  governs  the  jump  in  stress  across  shock  waves. 

The  Prony  series  approximation  to  the 
stress  relaxation  function  is  particularly  convenient 
for  numerical  implementation  in  explicit  codes,  as 
it  leads  to  a  simple  rate  form  for  the  shear  stress. 
For  a  given  material,  the  number  of  terms,  A,  that  is 
needed  in  (1)  depends  on  the  time  window  over 
which  the  viscoelastic  response  is  to  be  modeled. 
For  simulations  of  impact  and  shock  waves, 
relaxation  times  on  the  order  of  fractions  of  a  nano¬ 
second  may  be  needed.  On  the  other  hand, 
relaxation  times  more  than  a  decade  or  two  larger 
than  the  duration  of  the  simulation  will  yield 
exponential  terms  close  to  unity  for  all  times  of 
interest,  so  that  the  corresponding  terms  in  (1)  may 
be  lumped  into  the  equilibrium  modulus.  A  general 


rule  of  thumb  is  that  relaxation  times  should  be 
equally  spaced  on  a  log  time  scale,  with  at  least  one 
term  for  each  decade  of  time. 

The  viscoelastic  properties  of  the 
propellant  were  obtained  by  Dynamic  Mechanical 
Analysis  (DMA).  The  response  to  steady  state 
torsional  vibration  was  measured  at  six  different 
frequencies  from  0.1  Hz  to  30  Hz.  The  shear  strains 
were  on  the  order  of  0.1%  and  thus  well  within  the 
region  of  linear  viscoelastic  response.  Each  of  the 
fixed  frequency  tests  involved  a  slow  (about 
0.5°C/min)  temperature  sweep  from  -100°C  to 
25 °C,  spanning  the  glass  transition  temperature 
(roughly  -58°C).  The  tests  were  performed  by  Rob 
Jensen  of  ARL. 


FIGURE  3.  Storage  modulus  as  a  function  of 
temperature  at  six  fixed  frequencies. 


The  storage  modulus  G'  and  loss  modulus 
G"  are  obtained  from  the  components  of  stress  in 
phase  with  the  strain  and  out  of  phase  with  the 
strain,  respectively.  They  are  plotted  in  Figures  3 
and  4.  Observe  that  the  peak  in  the  loss  modulus  is 
an  order  of  magnitude  less  than  the  corresponding 
value  of  the  storage  modulus. 


temperature  for  six  fixed  frequencies. 


Another  advantage  of  the  Prony  series 
approximation  is  that  it  allows,  at  least  in  principle, 
the  straightforward  determination  of  the  moduli 
and  Gn  in  (a)  from  the  DMA  data,  since  in  this  case 
the  storage  modulus  G'  is  given  as  a  function  of 
the  radial  frequency  co  by  the  relation 

G'(«0  =  G.+iG.r^V  (3) 

tt  1  +  (jnO)) 


Note  that 

Gf(0)  =  G<o  and  G'(oo)  =  G(0)  (4) 

Assuming  that  the  number  of  terms  N  and 
corresponding  relaxation  times  Tn  has  been  chosen, 
one  may  determine  the  values  of  and  Gn  in  (3) 
that  best  fit  the  storage  modulus  data.  The  difficulty 
with  this  approach  is  that  the  DMA  data  does  not 
cover  a  wide  enough  frequency  range  to  accurately 
determine  the  moduli.  This  problem  is 

circumvented  by  using  time-temperature 
superposition  (TTS)  and  involves  several  steps. 

First,  the  fixed  frequency  data  in  Figures  3 
and  4  are  converted  to  modulus  vs.  frequency 
curves  at  fixed  temperatures.  Next,  a  reference 
temperature  T0  is  selected,  and  the  storage  modulus 
vs.  frequency  curves  at  the  various  temperatures  are 
shifted  horizontally  relative  to  the  curve 
corresponding  to  the  reference  temperature  so  as  to 
form  a  single  smooth  curve,  called  a  master  curve. 
The  principle  behind  this  procedure  is  the 
assumption  that  the  temperature  dependence  of  the 
stress  relaxation  function  G  is  governed  by  the 
simple  relation 

G{t,T)  =  G(t  /aT,T0)  (5) 

where  the  shift  factor  aT  is  a  decreasing  function  of 
the  temperature  T  with  the  value  1  at  T  =  T0  .  This 
implies  that  all  of  the  relaxation  times  scale 
proportionally  with  temperature: 

Zn(T)=ar  Tn(TQ)  (6) 

and  that  the  temperature  dependence  of  the  storage 
and  loss  moduli  are  given  by 

G\coJ)=G\aTco,T0)  , 


G"(co,T)=G"(aTco,T0)  (7) 

On  a  log  frequency  scale,  (7)  yields  the  relation 

G’(iog^,r)=G'(iog^  +  iogflr,r0)  (8) 

so  that  modulus  vs.  log  frequency  curves  at 
different  temperatures  superimpose  upon  shifting 
horizontally  by  a  temperature-dependent  factor. 
Here,  T0  was  chosen  to  be  the  glass  transition 


temperature,  -58°C.  The  resulting  storage  modulus 
master  curve  is  shown  in  Figure  5.  Note  that  the 
frequency  range  has  been  extended  from  about  2Vi 
decades  to  over  20  decades.  The  corresponding 
temperature-dependent  horizontal  shift  factors  are 
given  by  the  upper  curve  in  Figure  6.  These 
“optimal”  shift  factors  were  determined  so  as  to 
yield  the  smoothest  possible  curve. 


FIGURE  5.  Storage  and  loss  moduli  master 
curves  with  no  vertical  shifts,  using  the  optimal 
horizontal  shifts  for  the  storage  modulus. 


FIGURE  6.  Optimal  horizontal  shift  factors  for 
the  storage  and  loss  moduli  (no  vertical  shifts). 

If  the  principle  of  TTS  were  strictly  valid, 
as  reflected  in  the  relation  (5),  then  by  (7)  we  see 
that  the  same  horizontal  shift  factors  should  work 
for  both  the  storage  and  loss  moduli.  But  in  fact  the 
optimal  horizontal  shift  factors  for  the  loss  moduli 
are  given  by  the  lower  curve  in  Figure  6.  Since  the 
vertical  axis  is  logarithmic,  this  results  in  a 
maximum  difference  of  about  VA  decades  in 
frequency  (equivalently,  time)  between  these  shift 
factors.  When  the  optimal  storage  modulus  shift 
factors  are  used  to  generate  the  loss  modulus  master 
curve,  the  result  is  the  bell-shaped  curve  in  Figure 
5.  Note  the  large  scatter,  particularly  at  higher 
frequencies.  Likewise,  had  the  optimal  loss 


modulus  shift  factors  been  used  to  generate  the 
storage  modulus  master  curve,  there  would  have 
been  a  substantial  scatter. 

A  common  practice  is  to  overlook  these 
discrepancies  and  simply  fit  the  coefficients  G^  and 
Gn  in  (3)  to  the  smooth  storage  modulus  master 
curve,  using  the  corresponding  shift  factors  for  this 
curve  to  obtain  the  temperature  dependence  of  the 
relaxation  times  via  (6).  However,  Simon  and 
Ploehn12  argue  against  this.  In  fact,  they  conclude 
that  “master  curves  based  on  superposability  of  the 
loss  modulus,  rather  than  the  storage  modulus,  may 
lead  to  better  representations  of  true  response  of  a 
material.” 

One  possible  explanation  for  the 
discrepancies  between  the  storage  and  loss  modulus 
shift  factors,  is  temperature  dependence  of  the 
equilibrium  and  instantaneous  elastic  moduli  G^= 
G(oo)  and  G(0).  Observe  that  (5)  implies 

G(0, T)  =  G(0,T0),G(co,T)  =  G(co,T0)  (9) 

That  is,  the  instantaneous  and  equilibrium  elastic 
moduli  are  temperature  independent.  These 
conclusions  are  at  best  approximations.  It  turns  out 
that  even  small  corrections  for  the  temperature 
dependence  of  the  instantaneous  and  equilibrium 
moduli  can  lead  to  substantial  changes  in  the 
horizontal  shift  factors. 

The  correction  most  commonly  used  in  the 
polymer  literature  is 

G(cc,T)  =  yG(co,TQ)  ,  (10) 

where  the  equilibrium  elastic  modulus  is 
proportional  to  the  absolute  temperature.  This 
relation  is  supported  by  experiment  as  well  as  the 
kinetic  theory  of  rubber  elasticity.  By  (4)  we  see 
that  (10)  is  equivalent  to 

G'(0,T)  =  ^rG'(0,T0)  (11) 

Note  that  the  zero  here  is  co  =  0,  which  corresponds 
to  log  co  =  -oo.  Clearly,  when  (9)  holds,  no  amount 
of  horizontal  shifting  can  bring  the  storage  modulus 
vs.  frequency  curves  at  different  temperatures  to 
coincidence,  since  they  have  different  low 
frequency  asymptotes.  The  most  common 
procedure  used  to  correct  for  this  is  to  first 
vertically  shift  the  GX®,7)  and  G"  (co,7)  curves 
by  multiplying  them  by  the  factor  Tq/T. ,  which  in 
effect  cancels  out  the  temperature  dependence  of 
the  equilibrium  modulus.  Then  the  horizontal 
shifting  procedure  described  above  is  applied.  This 


generally  works  well  for  temperatures  above  the 
glass  transition.  When  applied  to  the  data  above, 
this  yielded  a  slight  improvement  in  the 
coincidence  of  the  storage  and  loss  shift  factors  at 
the  higher  temperatures,  but  resulted  in 
substantially  worse  discrepancies  at  the  lower 
temperatures.  This  is  due  to  the  fact  that  the 
temperature  dependence  of  the  instantaneous  elastic 
modulus  has  not  been  properly  accounted  for. 

In  contrast  with  the  equilibrium  modulus, 
the  instantaneous  elastic  modulus  G(0,7)  =  G'(°°,7) 
should  decrease  with  temperature.  McCrum  and 
Pogany13  and  Schapery14  have  emphasized  the  need 
to  account  for  this  when  performing  TTS  on  data 
below  the  glass  transition  temperature.  We  have 
essentially  followed  their  procedure  here.a  We  omit 
the  details,  but  the  basic  idea  is  as  follows.  First,  the 
appropriate  temperature  dependent  equilibrium 
modulus  as  determined  by  (10),  is  subtracted  from 
the  storage  modulus  vs.  frequency  curve  at 
temperature  T,  so  that  all  curves  have  a  low 
frequency  asymptote  of  zero.  However,  due  to  the 
dependence  of  G'(°°,7)  on  T,  the  high  frequency 
asymptotes  of  these  vertically  shifted  curves  will 
not  coincide,  so  a  second  multiplicative  scaling  is 
performed  to  correct  for  this.  We  assumed  that 
G'(oo,T)  varies  linearly  with  T.  Thus  the 
parameters  required  by  this  procedure  are 
G'(°°,Tb),  the  slope  p  =  dG\oo,T)/dT,  and  the  value 
of  G'(0,r0)  for  use  in  (1 1).  Initial  estimates  for  the 
latter  two  parameters  were  obtained  from  Figure  5. 
We  are  not  aware  of  any  data  from  which  the  slope 
p  could  be  inferred.  The  data  in  Figure  3  provide 
values  for  dG\<&,T)!dT  only  for  co  in  the  range 
from  0.1  to  30  Hz.,  whereas  P  is  the  limit  of  this 
derivative  as  co  — >  oo.  We  took  as  our  initial 
estimate  for  p  the  smallest  slope  on  the  modulus  vs. 
temperature  curves  in  Figure  3.  Vertical  shifts  of 
the  data  were  performed  as  described  above, 
followed  by  horizontal  shifts,  to  obtain  a  new  set  of 
storage  and  loss  modulus  master  curves. b  Several 
iterations  were  performed,  with  the  three 
parameters  modified  slightly  at  each  iteration.  The 
final  value  of  p  used  was  -0.004  GPa/°C,  which  is 


a  Their  corrections  are  for  the  case  where  the  data 
are  compliances,  but  an  analogous  procedure  works 
for  the  moduli. 

b  A  slightly  different  (in  fact,  simpler)  vertical 
shifting  procedure  is  required  for  the  loss  modulus. 


slightly  less  than  the  smallest  slope  of  the  curves  in 
Figure  3.  This  results  in  a  decrease  in  the 
instantaneous  elastic  modulus  of  about  0.5  GPa  as 
the  temperature  increases  from  -100°C  to  25°C. 

The  resulting  master  curves  and  shift 
factors  are  shown  in  Figures  7  and  8.  Observe  that 
the  optimal  horizontal  shifts  for  the  storage  and  loss 
moduli  are  now  nearly  identical  except  at  the 
highest  temperatures.  The  smooth  master  curves  in 
Figure  7  were  generated  from  the  same  shift 
factors,  which  are  essentially  the  average  of  those 
in  Figure  8.  The  coefficients  G^  and  Gn  in  (3)  were 
then  fit  to  this  storage  modulus  master  curve  using 
35  terms.  When  used  in  the  Prony  series  (1),  these 
yield  a  smooth  master  curve  for  the  stress 
relaxation  function  G(t,T0 )  that  covers  20  decades 
of  time.  A  portion  of  that  curve,  corrected  for  room 
temperature  using  the  appropriate  shift  factor  from 
Figure  8,  is  shown  in  Figure  9. 


FIGURE  7.  Storage  and  loss  moduli  master 
curves  with  vertical  shifts,  using  the  same 
optimal  averaged  horizontal  shifts. 

Also  shown  in  figure  7  is  an  8-term  fit  for  the  time 
window  spanning  0.25  ns  to  30  ms.  This  should  be 
more  than  sufficient  for  most  impact  simulations. 


FIGURE  8.  Optimal  horizontal  shift  factors  for 
the  storage  and  loss  moduli  after  vertical  shifts. 
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FIGURE  9.  Prony  series  approximations  to  the 
stress  relaxation  function  at  23°C. 

Damage  Parameters 

As  stated  earlier,  the  total  strain  rate  is  the 
sum  of  the  deviatoric  viscoelastic  strain  rate  and  the 
deviatoric  cracking  strain  rate.  This  generalization 
allows  ductile  or  brittle  material  behavior  in  the 
fracture  model.  The  failure  model  assumes  that  for 
each  element  a  penny-shaped  micro-crack  exists 
normal  to  the  direction  of  the  maximum  (principal) 
deformation  rate,  as  shown  in  Figure 
10. 


FIGURE  10.  The  viscoSCRAM  hot  spot  model 
showing  friction  generated  along  a  crack  face9. 


The  deviatoric  cracking  strain  rate  is  derived  as  a 
function  of  the  average  crack  radius,  c,  and  the 
initial  flaw  size,  a.  An  evolution  equation  for  the 
crack  growth  rate  is  assumed  to  depend  on  the 
stress  intensity  factor  as  shown  in  Figure  1 19. 


Stress  Intensity,  K=(7tcg)0'5 


FIGURE  11.  Schematic  of  crack  growth  rate,  c . 


A  rate-dependent  damage  function,  vmax,  is 
used  to  calibrate  the  accumulation  of  internal 
damage  and  appears  to  be  the  most  important 
parameter  in  the  damage  model.  The  threshold  for 
stress  intensity,  K0,  is  a  function  of  the  coefficient 
of  static  friction,  pst,  and  a  stress  intensity 
parameter,  m.  K0  is  also  related  to  the  effective 
deviatoric  strain  rate  (at  each  time  increment)  in 
order  to  convert  from  a  tensor  form  to  an  equivalent 
unidirectional  strain  rate.  This  is  why  the  current 
version  of  viscoSCRAM  is  isotropic.  The  values  for 
K0,  m,  and  pst  are  varied  from  material  to  material. 
Since  vmax  controls  the  yield  strength  of  the 
material,  it  is  necessary  to  optimize  vmax  to  at  least 
three  different  compressive  or  tensile  loading  rates. 

Hot  Spot  Ignition  Model 

Thermal  heating  in  viscoSCRAM  includes 
bulk  heating  at  the  continuum  length  scale  and  hot 
spot  heating  at  the  microstructural  length  scale. 
Bulk  heating  includes  mechanical  terms  describing 
viscous,  damage,  and  adiabatic  volume  change  as 
well  as,  a  chemical  decomposition  term.  Chemical 
decomposition  is  based  on  Arrhenius  first  order 
chemical  kinetics.  For  the  continuum,  the  rate  of 
temperature  change  with  respect  to  time  is  written 

t  =  -  rTsu +—[(«').,  +M,  ]+fU,  d2) 

P°v 

where  the  first  term  on  the  right  hand  side 
represents  adiabatic  compression  heating  rate,  the 
second  term  represents  the  inelastic  work  rates  due 
to  viscoelastic  effects  and  cracking  damage,  and  the 
third  term  is  the  bulk  chemical  heating  rate. 

The  ignition  criterion  in  the  viscoSCRAM 
hot  spot  model  describes  frictional  heating  due  to 
crack  faces  sliding.  Given  the  stresses  from  two 
adjacent  elements,  the  local  strain  energy  release 
rate  is  determined.  Then,  at  the  end  of  a  time  step  in 
the  simulation,  the  change  in  crack  length  of  the 
interface  crack  is  determined.  If  the  interface  crack 
grows  to  be  wider  than  the  length  of  the  element 
edge,  the  interface  fails  and  is  allowed  to  separate 
by  not  enforcing  the  constraints  on  the  adjacent 
interface  nodes.  As  the  simulation  progresses,  the 
failed  interfaces  coalesce  into  macroscopic  cracks. 
Once  the  shear  stress  exceeds  a  slip  criterion,  the 
adjacent  crack  faces  are  assumed  to  slip.  The  work 
done  by  the  slipping  faces  will  generate  heat  and 
possibly  ignite  the  energetic.  This  frictionally 


triggered  hot  spot  model  is  included  in  the  energy 
balance  on  a  differential  material  volume  near  the 
crack  face  along  with  mechanical  and  chemical 
heating  terms  (contribution  of  heat  conduction  can 
be  neglected).  Referring  to  Figure  10,  the  heat 
transfer  near  the  1  -D  crack  face  is  given  by 

_E^ 

PfCft=kfTu+pfAHZe  RT+pstpeij,  lf  >  y  >  0  (13) 

_E^ 

pscsT=ksTii+psAHZe~™,  y>lf  (14) 

where  the  hot  spot  length  scale  is  lj  and  just  is  the 

coefficient  of  static  friction.  In  eqns.  (13)  and  (14), 
the  left  hand  side  is  the  heat  stored  in  the  region  of 
the  hot  spot.  The  first  term  on  the  right  side  is  the 
heat  conducted  away  from  the  hot  spot  and  the 
second  term  is  the  chemical  heat  generation  per  unit 
volume.  For  each  finite  element,  the  deviatoric 
stress  is  found  on  a  plane  normal  to  the  direction  of 
the  maximum  principal  deformation  rate.  If  the 
maximum  shear  stress  exceeds  the  value  of  justp 
then  the  crack  is  assumed  to  slip  and  generate  heat. 
Note  that  p  is  the  compressive  pressure  and  if  it  is 
positive  the  crack  is  open  and  will  not  generate 
heat. 

COMPUTATIONAL  RESULTS 

The  simulations  of  the  shear  punch 
validation  test  were  performed  on  a  parallel 
computer  platform  at  the  ARL  MSRC  using  the 
LLNL  code  ALE3D15.  The  entire  computational 
domain  included  the  incident  and  output  bars,  the 
specimen,  and  the  specimen  holder.  For  the 
simulations  presented,  the  50cm  striker  bar  was 
replaced  with  a  prescribed  input  velocity  boundary 
condition  on  the  end  nodes  of  the  incident  bar.  The 
z- velocity  pulse  had  a  1300m/s  material  velocity,  a 
5ps  rise  time  with  duration  of  200 ps. 

A  hybrid  computational  domain  was  also 
built  for  the  simulations  using  8-node  hexagonal 
elements.  Slide  surfaces  and  symmetry  conditions 
were  also  used  to  create  the  ^-symmetry,  butterfly 
computational  domain,  as  shown  in  Figure  12.  The 
input  bar,  output  bar  and  specimen  holder  were 
modeled  using  an  elastic-plastic  description.  The 
specimen  constitutive  behavior  was  described  using 
the  viscoSCRAM  model. 


FIGURE  12.  The  hybrid  computational  domain 
used  for  the  shear  punch  simulations. 

A  plot  of  the  deformed  mesh  for  the 
specimen  is  shown  in  Figure  13.  Note  the  localized 
shear  surface  that  formed  in  the  specimen.  The 
localized  strain  in  the  specimen  emanated  from  the 
periphery  of  the  indenting  piece  and,  later  in  time, 
formed  on  the  distal  end  at  the  holder/output  bar 
interface.  Figure  13  also  shows  a  plot  of  the 
surrogate  specimen  temperature  at  the  end  of  the 
simulation.  PC  has  a  melt  temperature  of  558  K. 
The  temperature  rise  is  due  to  the  conversion  of 
plastic  work  to  heat.  Although  the  temperature 
localizes  near  the  bar/specimen  interface,  it 
dissipates  to  neighboring  elements  because  of  the 
mesh  resolution.  For  a  finer  mesh,  the  temperature 
may  localize  along  the  idealized  shear  surface  and 
reach  a  higher  order  of  magnitude. 

The  specimen  geometry  is  different  from 
what  is  required  in  a  conventional  Kolsky  bar.  For 
this  reason  the  strain  rate  is  not  uniform  in  the  shear 
punch  specimen.  The  specimen’s  strain  rate  reaches 
~8000-9000s_1  and  localizes  along  the  idealized 
shear  surface.  An  examination  of  the  shear  stress  in 
the  specimen  during  compressive  loading  at  600ps, 
shows  the  stresses  reach  40-50  MPa.  By 
comparison  the  principal  compressive  stress 
reaches  ~150MPa  in  the  center  region  and 
~300MPa  in  the  outer  region.  Of  course,  the  state  of 
stress  in  the  specimen  will  change  at  the  arrival  of 
the  transmitted  wave.  For  a  finer  specimen  mesh 
resolution  subjected  to  this  complex  state  of  stress, 
the  specimen  material  may  form  adiabatic  shear 
bands.  We  also  note  that  the  pressure  in  the  PC 
specimen  reaches  approximately  -5MPa  (tensile 
hydrostatic  stress).  This  pressure  is  above  the 
fracture  pressure  (-80MPa)  therefore,  the  PC 
specimen  did  not  fracture  in  this  simulation. 
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FIGURE  13.  Final  deformed  shape  and 
temperature  [K]  in  the  surrogate  PC  specimen. 

A  comparison  of  the  strain  gage  signals  to 
the  observed  result  shows  excellent  agreement,  as 
shown  in  Figure  14.  The  incident  and  reflected 
pulses  are  shown  in  Figure  14a.  The  curvature  at 
the  beginning  of  the  experimental  input  pulse  is  due 
to  wave  shapers  added  in  front  of  the  input  bar. 
There  were  no  wave  shapers  added  in  the  numerical 
simulations.  The  ringing  seen  at  the  beginning  and 
end  of  the  numerical  incident  signal  are  due  to  the 
sharp  discontinuity  of  the  prescribed  velocity 
boundary  condition.  Smoothing  this  boundary 
condition  will  reduce  the  ringing. 
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FIGURE  14.  Comparison  of  strain  signals  for 
the  incident  bar  strain  gage  and  the  transmitted 
strain  signal. 


The  viscoSCRAM  constitutive  model 
described  was  used  to  represent  the  behavior  of  the 
double-base  propellant  specimen  in  the  dynamic 
shear  punch  test,  as  shown  in  Figure  15.  New  crack 
faces  are  created  during  the  early  loading  stages.  As 
a  result,  Figure  15  illustrates  that  the  propellant 
generates  heat  due  to  chemical  decomposition 
shortly  after  the  arrival  of  the  dynamic  compression 
wave.  In  the  experiment,  the  double  base  propellant 
experienced  more  plastic  deformation  and  cracking 
(resembled  an  extrusion  process)  before  it 
generated  heat  from  chemical  decomposition. 


FIGURE  15.  Evolution  of  chemical 
generation  using  the  viscoSCRAM  model. 
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SUMMARY  AND  CONCLUSIONS 


Numerical  simulations  of  a  shear  punch 
test  have  been  completed  to  study  the  effects  of 
shear  loading  on  various  energetics.  To  date  we 
have  completed  simulations  for  nonenergetic 
polymer  materials,  plastic  bonded  explosives,  and 
double  base  propellants.  For  viscoelastic 
parameters,  vertical  and  horizontal  shifts  of  DMA 
data  should  be  performed  to  obtain  a  set  of  storage 
and  loss  modulus  master  curves.  For  the  damage 
model,  the  parameters  that  control  the  yield 
strength  must  be  compared  and  fitted  with  observed 
data  to  obtain  a  damage  law.  Simulations  showed 
excellent  agreement  of  the  strain  gage  signals  and 
showed  the  general  trend  of  an  idealized  shear 
surface  in  the  specimen.  The  hybrid  mesh 
capability  enabled  complete  modeling  of  the  shear 
punch  test.  The  Lagrangian  formulation  used  for 
the  incident  bar  and  output  bar  provided  an  efficient 
solution  to  wave  propagation.  The  ALE  mesh  for 
the  specimen  prevented  hourglassing  and  excessive 
material  advection  while  maintaining  a  reasonable 
timestep.  More  work  is  needed  to  reduce  the 
advection  in  the  specimen  for  the  simulations,  i.e. 
make  the  specimen  more  Lagrangian. 

The  hot  spot  shear  initiation  model  was 
included  in  viscoSCRAM  for  a  double-base 
propellant.  The  simulation  predicted  chemical  heat 
generation  at  the  early  stages  during  the  arrival  of 


the  dynamic  compression  wave.  Further  work  is 
required  on  determining  the  sensitivity  of 
viscoSCRAM  input  parameters.  Furthermore,  the 
authors  believe  that  the  isotropic  behavior  of  the 
damage  function  causes  cracks  to  accumulate  in  all 
directions.  Zuo  et  al.16  recently  showed  the 
importance  of  stable  and  unstable  orientation  on 
shear  cracks.  Since  the  viscoSCRAM  model  is  a 
“work  in  progress”,  a  re-formulation  of  the  damage 
function  into  a  tensor  quantity  would  directly 
influence  ignition  criteria.  This  would  progress 
viscoSCRAM  into  a  very  useful  model  for 
predicting  insensitive  munition  behavior. 
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